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Multigrid  Applied  to  Singular  Perturbation  Problems 


David  Kamowitz* 


Abstract 

The  solution  of  the  singular  perturbation  problem 

—  eu"  4-  b[x)u'  =  /,  0  <  i  <  1 


1  3>  £  >  0,  u(0)  =  u0l  u(l)  =  uj 

by  a  multigrid  algorithm  is  considered.  Theoretical  and  experimental  results  for  a  number  of 

different  discretizations  are  presented.  The  theoretical  and  observed  rates  agree  with  the  results 

developed  in  an  earlier  work  of  Kamowitz  and  Parter. 

In  addition,  the  rate  of  convergence  of  the  algorithm  when  the  coarse  grid  operator  is  the 

natural  finite  difference  analogue  of  the  fine  grid  operator  is  presented.  This  is  in  contrast  to 

the  case  in  the  previous  work  where  the  Galerkin  choice  (iff  L^lfa)  was  used  for  the  coarse  grid 

operators.  - - 
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1  Introduction 


This  report  discusses  the  application  of  a  multigrid  algorithm  to  the  solution  of  the  following  one 
dimensional  singular  perturbation  problem: 

—  eu"  +  b(x)u'  =  /,  0  <  x  <  1  (l) 

with 

1  »  e  >  0,  u{0)  uo,  u(l)  =  uj. 

Many  other  authors  have  discussed  the  application  of  various  methods  of  solution  to  the  alge¬ 
braic  problem;  in  particular  see  Dorr  [Dorr70a],  Babuska  [Babu69a],  Ervin  and  Layton  [Ervi85a], 
and  Kellogg  and  Tsan  [Kell78a]. 

Much  of  the  literature  regarding  multigrid  methods  restricts  itself  to  the  solution  of  nice  prob¬ 
lems.  Indeed  most  authors  require  that  the  linear  system  be  well  conditioned  in  addition  to  symmet¬ 
ric  and  positive  definite.  However,  in  the  case  of  singular  perturbation  problems  as  the  coefficient  of 
the  second  order  term  tends  to  zero  the  usual  symmetric  discretization  fails  to  be  of  positive  type. 
Thus  the  first  measure  taken  in  the  numerical  discretization  of  these  problems  is  to  replace  the  usual 
symmetric  difference  of  the  first  order  term  with  some  form  of  skewed  differencing.  In  particular 
for  problems  with  turning  points  this  may  lead  to  a  system  of  equations  which  is  ill-conditioned 
for  small  e. 

From  the  standpoint  of  calculating  a  numerical  approximation  to  the  solution  of  problem  (1) 
the  first  question  is:  does  the  discretization  converge  to  the  continuous  solution?  Then,  assuming 
it  does,  how  does  the  multigrid  algorithm  perform  as  a  solver  for  this  system  of  linear  equations? 
What  modifications,  if  any,  are  necessary  in  the  multigrid  algorithm? 

The  main  result  shows  that  if  the  original  system  of  equations  resulting  from  the  discretization 
of  problem  (1)  is  of  positive  type  then  the  theoretical  results  for  the  multigrid  algorithm  developed 
in  Kamowitz  and  Parter  [Kamo85a]  and  in  McCormick  and  Huge  [McCo82a]  apply. 

From  a  computational  perspective  it  is  convenient  to  use  for  the  coarse  grid  operators  of  the 
multigrid  algorithm  the  operators  that  are  the  finite  difference  analogue  of  the  original  operator. 
In  Section  6  the  rate  of  convergence  of  the  algorithm  using  the  finite  difference  version  of  the  coarse 
grid  operators  is  considered.  It  is  shown  that  the  new  rale  of  convergence  is  an  0(h 2)  perturbation 
of  the  rate  obtained  using  the  Galerkin  choice. 


1 


2  The  Discrete  Problem 


Three  model  problems  are  chosen  for  study,  one  where  the  sign  of  6(x)  does  not  change  and  two 
where  b[x)  changes  sign.  The  three  problems  are  designated 


Problem  BL 


Lu  —  —eu"  +  u1  —  0,  0<x<l 


Problem  TP-1 


Lu  =  -  eu"  +  (x  -  —  )u'  =  0,  0  <  x  <  1 

2 


Problem  TP-2 


Lu  =  —eu"  —  (x  —  —  )u;  =  0  0  <  x  <  1. 

v  2 


(2) 

(3) 

(4) 


For  all  three  problems  the  boundary  conditions  are 

tt(0)  =  l,  «(1)  =  3. 

For  the  discrete  problem  as  usual  let  N  >  0  be  chosen  and  set  h  —  1  /(N  +  1).  The  interval 

n  =  (o,i) 

is  discretized  to  form 


Qh  =  {ih  :  1  <  t  <  N}. 


The  notation  x,  refers  to  the  point  ih  and  u,  refers  to  u(ii).  In  addition  the  usual  notation  for 
finite  differences  is  used: 


ti,  +  l  -  u,  U.-U.-j 

D+U,  =  - : - ,  D-Ui  =  - r - 


n  «,+  I  “«.-l  ,,  n  «i+l  -  2«i  +  «,•_! 

IT)U,  -  - — - — ,  D+D-Ui  - 


2  h  ’  - -  h2 

For  problem  ( BL)  the  following  two  discretizations  are  considered: 


L\ui  r:  -eD+D_Uj  +  D_u,  (upwind  differencing) 


— D+D-ut  f  D_u,  (see  Kellogg  and  Tsan). 


Since  problems  (TP-1)  and  (TP-2)  have  turning  points  at  x  =  j  it  is  necessary  to  change  the 
direction  of  the  discretization  of  the  first  order  term.  For  problem  (TP-1)  the  discretization  tested 
is 

L\\n  = -eD+D-Ui-k- (xi  -  \)D+Ui,  1  <  i  <  (5) 

L\ui  =  - eD+D-Ui  +  (xi  -  +  *  +  1  <  «  <  N.  (6) 

Similarly,  for  problem  (TP- 2)  the  discretization  is 

L\ui= -eD+D-Ui- (ii-  \)D-Ui,  1  <  t  <  (7) 

L\v.i=  -eD+D-Ui- (xt-  \)D+  uit  — +  1  <  i  <  N.  (8) 

Note  that  each  of  the  discretizations  L*,  k  =  1 ,  ...  ,4  results  in  a  tridiagonal  linear  system  of 
equations  whose  coefficient  matrix  may  be  denoted  by 

Lkh  =  [  -<*i  A  ~Tx  }, 

with 

a,  >  0,  7,-  >  0,  and  /?,  >  a,  +  7.- 

Thus  each  L k  —  1, ...  ,4  is  an  M-matrix  and  the  linear  system  of  equations 

Lkhui  =  /,,  1  <  t  <  N, 

u o  and  tijv+i  fixed  has  a  solution;  see  for  example  Berman  and  Plemmons  [Berm79a], 

3  Review  of  the  Multigrid  Algorithm 

The  particular  rnultigrid  algorithm  used  to  solve  (1)  has  been  discussed  in  detail  in  Kamowitz  and 
Farter  [Kamo85a];  thus  only  a  cursory  review  is  given  here.  In  particular  the  details  of  the  theory 
behind  the  convergence  results  can  be  found  there.  What  is  important  to  realize  is  that  although 
the  algorithm  and  convergence  results  discussed  in  the  previous  work  apply  to  well  conditioned  two 
point  boundary  value  problems  (he  same  bounds  on  the  rate  of  convergence  apply  to  the  singular 
perturbation  problems  discussed  here. 

In  order  to  completely  describe  the  algorithm  a  number  of  spaces  and  operators  need  to  be 


where 


p  =  29 

Denote  by  Skh  the  space  of  grid  functions  defined  on  Clkh-  From  now  on  the  notation  kh  will  be 
replaced  by  k.  Denote  by  Gk  the  smoother,  /£+J  the  restriction  map,  and  the  interpolation 

map.  Associated  with  each  space  Sk  there  is  also  a  nonsingular  tridiagonal  operator 


Lk  :  St  ►->  Sk, 


where  Lk  is  of  positive  type.  The  operators  Lk  are  again  denoted 

Lk  =  |  P*  -7?  ]• 

These  operators  wiil  be  formally  defined  after  the  statement  of  the  algorithm. 

The  following  is  an  outline  of  the  algorithm  used.  Assuming  t/”  is  the  nth  approximation  to 
the  solution  of  the  system  of  equations 

LxUx  =  Fi 

algorithm  MG(U£,  Lk,Fk,k)  returns  U£+1,  the  next  iteration  of  the  multigrid  algorithm.  The  grid 
layer  is  denoted  by  k\  set  A:  =  1  to  start  the  algorithm  on  the  finest  grid  layer. 


Algorithm  MG(U£,  Lk,  Fk,  k) 

1.  Coarse  Solve:  If  k  =  g  (coarsest  grid)  then  return 

MG*-  Lg-lFg 

otherwise 

2.  Smooth:  Apply  the  smoother,  Gk,  call  the  result  of  this  step  U£. 

If  k  —  1  use  U['  as  the  initial  guess,  otherwise  use  0. 

3.  Recursively  Apply  the  Algorithm:  Set 


U 


n  *-  1 


0?  +  j£+lA/G(0,  Lk+i,  lV  V*  -  LkU?),  k  +  1). 


4,  (optional)  Smooth  Again:  Set  Un+I 

5.  Return:  Set  MG  *-  f/n+1. 


GkU 


k  . 

nil 


As  defined  algorithm  MG  is  called  a  symmetric  V-cycle  if  step  4  is  used;  otherwise  algorithm 
MG  is  referred  to  as  a  slash  cycle  (following  the  notation  of  McCormick  and  Ruge). 

For  the  smoother,  G,  choose  m-applications  of  damped  Jacobi  iteration  with  parameter  a 
Formally,  repeat  for  1  <  r  <'  rn 

I  he  interpolation  operator.  1(  is  defined  as  follows.  For  points  common  to  Uk  f  i  and  to  f 


and  at  odd  (new)  points  of  fi*  require 


{M4\i£/]}2j-i  =  °- 


This  results  in  the  explicit  system  of  equations  at  the  point  X2j-i 


Vti-i  =  J ■ 

Plj-  1 


The  restriction  operator,  /£fl,  is 


[Ikk  +  lU% 


a 


2;  77* 


-yfc 
72j 


TTK  a-  77*  4-  '■‘-1  rr* 

it - U2>-1  +  U2j  +  -ok  U2j+l 


@2j+l 


Note  that  if  Lk  is  symmetric  then 


A  fundamental  observation  due  to  McCormick  and  Ruge  [McCo82a]  in  the  symmetric  case  is  that 
S*  can  be  written 

Sk  —  Range  Jk+1  ©  Nullspace  /£+1  Lk.  (9) 

For  non-symmetric  problems  the  above  decomposition  follows  directly  from  the  characterization  of 
Range  /£+1  and  Nullspace  Ik+1Lk ;  see  Kamowitz  and  Parter. 

Finally  the  coarse  grid  operators,  Lk+i,  are  chosen  by  setting 


I'k+i  —  Lk+ 1  =  Ik+1  LkIk+l. 


A  direct  computation  shows  that 


t  _  f  fc  +  l  nk~ fl 

k'k  +  l  —  (  ~<*j  Pj 


„*+l 


with 


a 


it+i 


a2ja2;-l 


[  Ki-X 


k+  1 


nk  a&2j-l  ^;Q2,+  1 
‘?2;  -  -  -  " 


@2)-  1  ^2j  (  1 


ict  1 


-  k 

9r2j9,2y+  l 


^2?  t 


2]  f  1 


(10) 

(11) 

(12) 


Note  that  the  choice  of  /,*.  i  is  the  ‘natural’  choice;  however  it  is  not  the  only  feasible  choice  for 


3.1  Convergence  Theory 


For  completeness  some  results  from  [Kamo85a]  are  repeated. 


f"  -  Utr ue  ~  Un 

be  the  error  in  the  nth  iterate.  Here  Utruc  corresponds  to  the  true  solution  of  the  algebraic  problem 

Uh  Utruc  —  /• 

The  rate  of  improvement  of  the  nth  iterate  of  algorithm  MG  is  then 


To  estimate  p ,  the  asymptotic  rate  as  n  — >  oo,  one  needs  to  bound  the  pn . 


For  later  use  define 


\x\\lh  =  (Lhx,x)  =  X]  [Lhx)jXj. 


First  the  two  grid  process  is  considered.  Given  an  initial  error 


C°  -  Utrue  ~  U°, 


the  two  grid  process  yields: 


1.  Smoothing 


c°  — ►  e°  —  G'c0 


where  G1  corresponds  to  the  linear  part  of  the  smoother  G.  Note  that  from  (9) 


where 


(°  -  T]  +  I*hU)2h 


r]  G  Nullspace  llh Lf,  and  w2h  G  S2h. 


2.  Restricting  the  residual  r ^  —  L/,?0  to  yields 


n2h  =  llhLhl°  =  llhLhTi  +  I2hhLhI»hw2h  =  1^  LhI2hw2h 


t]  G  TV ullspacc  l\h  A/,. 


m 


-r.v; 


3.  Computing  the  coarse  grid  correction  directly  (as  is  the  case  for  the  two  grid  algorithm)  is 
equivalent  to  solving 

Llh'l’lh  =  Rlh  =  L2h^2h- 

Here 

Lih  =  Uh  =  l\hLkI$h 

SO 

*p2h  =  U’2A- 

4.  Finally,  correcting  c°  using  the  coarse  grid  correction  yields  e1: 

f1  =  Utruc-Ul 

=  Utrue  ~  {0°  +  Iihfrh) 

=  e°-$hihh 

=  T]  +  ^2hw2h  ~  I-lh*p2h 

-  T,. 


Note  that  if 


f°  €  S2h 


then  one  step  of  the  two  grid  procedure  solves  the  problem!  In  general,  since  ||G'||  <  1  the  rate  pl 
satisfies 


P  = 


< 


In  this  case 


so 


Nl 


Notice  that  the  rj  term  is  related  to  the  action  of  the  smoother  while  the  I^h w^h  term  is  eliminated 
by  the  multigrid  process  itself. 

In  Kamowitz  and  Farter  [Kamo85aJ  an  explicit  decomposition  of  was  found  in  terms  of  the 
eigenvalues  and  eigenvectors  of  the  damped  Jacobi  scheme.  This  decomposition  was  exploited  to 
compute  bounds  on  the  rate  of  convergence. 

For  the  two  grid  scheme,  the  bounds  on  p,  for  a  given  rn  and  a,  are  given  in  I  able  1.  Note  that 


a 

.333 

.500 

.667 

.750 

1.00 

1.333 

.500 

.333 

.400 

.429 

.500 

.572 

.250 

111 

.160 

.184 

.250 

.326 

.125 

.078 

.088 

.093 

.125 

.187 

.068 

.062 

.068 

.072 

.083 

.109 

Table  1:  Predicted  Rates  for  2  Grids 


.333 

.500 

.667 

.750 

1.00 

1.333 

1 

.633 

.577 

.561 

.561 

.577 

.614 

2 

.435 

.408 

.417 

.424 

.447 

.475 

3 

.336 

.335 

.349 

.357 

.378 

.403 

4 

.283 

.293 

.307 

.314 

.333 

.357 

Table  2:  Predicted  Asymptotic  Rate 

the  optimal  rate  is  obtained  for  a  —  .5.  In  the  succeeding  sections  the  bound  in  the  a  —  .5  case  will 
be  obtained  experimentally  for  problems  (BL),  (TP-1),  and  (TP-2). 

In  the  case  where  the  number  of  grids  is  arbitrary,  estimates  based  on  the  ideas  in  [McCo82a] 
are  given  in  Table  2.  As  indicated  in  jKamo85a|  these  bounds  are  not  sharp. 

4  The  Singular  Perturbation  Case 

For  the  study  of  singular  perturbation  problems  it  is  necessary  that  if  L*  is  an  M-mn'  rix  then  the 
coarse  grid  operators  L* ,  i ,  k  1,2,.  . .  are  also  M-matrices.  This  is  necessary  to  insure  that  each 
subproblem 

Bk<r\4>kY\  =  fk+  1 

is  solvable.  In  particular  since  the  smoother  Gk+  i  depends  on  Lk+\<  then  if  Lk+  i  is  an  M -matrix, 

!|g*+i||  <  i 


Lemma  4.1  If 


is  an  M-Matrix  then 


Ik  ••  [  of  p)  7*] 


r  I  ~  + 1  ok  f 1  k  -f  1  1 

I'kY  1  !  -<*j  Pj  -  lj  I 

is  also  an  M -matrix,  where  o^41,  /?*"*  ’  and  7* 4  1  arc  given  by  equations  (10  12). 


I 


1 


v  •  *  v 


V 


i''  A 


8 


Proof:  The  hypotheses  that  L *.  is  an  M-matrix  insures  that 


a*  >0,  /3k  >  0,  and  7*  >  0 


so  a*+1  and  7^+1  are  also  positive.  In  addition,  /3k+l  must  satisfy 

/?■  rl  >  «-  +  l  +  7*+1- 


Recalling  the  expressions  for  ak,  and  7^  note  that  this  is  equivalent  to  requiring 


^•Q2;  +  l  >  <*2ja2j-l 
02j+\  @2j—  1  @2j+ 1 


or, 


a" 


2j 


a 


2>-l 


i2]~ 


+  [a2;  +  ]  +  72,+ 1] 


2j+  1 


which  follows  directly  from  the  fact,  that  L *  is  diagonally  dominant. 

Since  each  of  the  Lk\-\  is  an  M-matrix,  the  theory  developed  in  Section  3.1  can  be  applied  in 
the  singular  perturbation  case. 


5  Experimental  Results 

The  algorithm  of  the  previous  sections  was  applied  to  problems  ( BL ),  (TP-1)  and  (TP-2).  Prob¬ 
lem  (PL)  exhibits  a  boundary  layer  at  1  -  1.  From  a  computational  point  of  view  the  system  of 
linear  equations  that  is  solved  is  well  conditioned  even  for  small  e;  however  the  fact  that  the  linear 
system  is  not  symmetric  leads  to  computational  difficulties  in  computing  the  experimental  rate  of 
convergence. 

In  the  case  of  problem  (TP- 1)  t  here  are  two  boundary  layers;  at.  x  ~  0  and  at  x  —  1 .  In  addition 
the  system  of  equations  that  is  solved  is  ill  conditioned. 

The  discrete  problem  corresponding  to  problem  (TP-2)  is  well  conditioned.  The  fact  that  there 
is  an  interior  turning  [mint  at  s  1/2  does  not  appear  to  lead  to  computational  difficulties. 

By  using  the  one  sided  diflerenee  schemes  discussed  in  the  previous  section  the  linear  systems 
arising  from  the  discretization  of  the  problems  are  of  positive  type. 

From  a  practical  standpoint,  however,  what  ellect  does  the  ill  conditioning  have  on  the  observed 
solution?  In  particular,  how  should  the  rate  of  convergence  (reduction  in  the  error)  be  measured? 


5.1  General  Remarks  Regarding  the  Experiments 

In  all  the  experiments  N ,  the  number  of  points  on  the  fine  grid,  was  63.  The  initial  guess,  f/°,  was 
constructed  so  that  the  initial  error  could  be  chosen  in  advance.  In  other  words, 

U°  =  u(Xj )  -  (y 

For  the  experiments  discussed  here 

The  number  of  smoothing  steps  for  the  experiments  discussed  here  was  set  to  one.  Experi¬ 
ments  were  run  where  the  number  of  smoothing  steps  was  greater  than  one.  When  more  than  one 
smoothing  step  was  used  there  was  no  observed  qualitative  difference  in  the  behavior  of  the  algo¬ 
rithm  compared  to  the  behavior  for  non-singular  perturbation  problems.  Unless  otherwise  noted 
the  results  are  for  the  two  grid  case.  The  computer  used  was  a  VAX-1 1  /780  and  the  tests  were  run 
in  double  precision  (roughly  16  decimal  digits  of  accuracy). 

The  damped  Jacobi  parameter,  a,  was  chosen  to  be  .5.  The  following  heuristic  argument  due 
to  Brandt  [Bran77a]  suggests  why  a  =  .5  is  optimal  for  one  smoothing  step.  His  suggestion  is 
to  choose  a  so  that  the  range  of  eigenvalues  which  are  reduced  by  the  damped  Jacobi  process  is 
equal  to  the  range  which  is  left  alone  by  the  process.  As  noted  in  [Kamo85a]  the  eigenvalues  of  the 
damped  Jacobi  scheme  come  in  pairs  A(p)  and  A (ft)  where 

u  \  V +  a  i  w  \  a  ~  (l 
AO,).  —  and  A(„)=  — 

and  ft  are  the  eigenvalues  of  the  scheme  for  a  -  0.  The  eigenvalues  //  are  real  and  distinct  and  as 
h  0  they  fill  out  the  interval  [0,  1],  Brandt’s  requirement  corresponds  to  choosing  a  so  that 

A ( -  I )  -  A(0) 

or  in  other  words  to  taking  a  1  /‘l.  Note  that  this  is  equivalent  to  requiring 

A(l)  A(0). 

Indeed  in  [Kamo8.r>a]  for  1  smoothing  step  the  theoretical  and  experimental  results  indicate  that 
this  choice  of  a  is  optimal. 

5.2  Results  for  the  Boundary  Layer  Problem 

For  problem  (HI,)  the  solution  to  the  analytic  problem  can  be  calculated  explicitly.  For  all  e  >  0 
t  he  solut  i< vn  <i,(.r)  is 
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Table  3:  Condition  Number  -  Problem  BL 


The  constants  C\  and  C2  satisfy 


Cx 


4/«  _  1  ’ 


C2  =  1  -  cx 


°'!c  -  3 


4/e  _  j  ' 


Note  that 


lim  Ci  =  0,  lim  C2  —  1 
6— »0  £— »0 


SO 


lim  ue(x)  =  1. 


e— *0 


Denote  by  f7/,t(x)  the  numerical  solution  of  problem  (BL)  for  a  fixed  h  and  e.  The  consistency 
condition  on  the  solution  then  requires  that  for  fixed  e 


lim  Uh  €(x)  =  ue(x). 

h — >0 

Dorr  [Dorr70a]  proved  that  for  fixed  h 

lim  Uh'C(x)  -  1. 

The  discretization  of  problem  (BL)  using  Llh  (standard  upwind  differencing)  results  in  an  ap¬ 
proximation  to  the  true  solution  which  is  an  0(h)  approximation.  To  improve  on  this  approximation 
Kellogg  and  Tsan  point  out  that  the  use  of  L\  results  in  an  0(h2)  approximation  for  e  >  0.  For 
the  reduced  problem  (e  —  0)  L2h  gives  an  O(h)  approximation. 

It  is  important  to  note  that  the  linear  system  of  equations  arising  from  the  discretization  of 
Problem  (BL)  is  not  ill-conditioned.  Table  3  displays  the  LINPACK  [Dong79aj  estimate  RCOND 
(an  estimate  of  the  inverse  of  the  condition  number)  for  L\,  h  -  1/64  and  e  =  1,  .01,  .005,  .003, 
and  .001. 

First,  some  general  conclusions  about  the  experimental  results.  For  both  L\  and  L\  the  al¬ 
gorithm  converged  to  the  solution  of  the  algebraic  problem  with  an  asymptotic  convergence  rate 
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identical  to  the  rate  predicted  in  Section  3.1.  Moreover  there  was  no  observed  qualitative  difference 
in  the  behavior  of  the  algorithm  with  respect  to  the  choice  of  L\  or  L\.  However,  as  e  — *  0  the 
behavior  of  the  iterates  changes  dramatically. 

Define 

IrJI 


*n  = 


rn-  1  [ 


where 

r„=  F  -  LhUn 

is  the  residual  after  the  nth  iteration  and  the  norm  used  is  the  l2  norm.  From  a  computational 
point  of  view  this  is  a  convenient  measure  of  the  rate  of  improvement  since  one  does  not  know  the 
true  solution.  From  past  computational  experience  this  ratio  is  bounded  above  by  the  theoretical 
error  reduction  rate.  Indeed  for  large  e,  say  e  =  1,  this  is  the  case.  For  small  e,  say  e  —  .001,  after 
a  small  number  of  iterations  Rn  exceeded  the  rate  predicted  in  Table  1,  and  then  as  n  increased  Rn 
declined  towards  the  predicted  rate;  see  the  solid  line  in  Figure  1.  The  dashed  line  will  be  referred 
to  later.  In  both  cases  two  grids  are  used. 

An  explanation  for  this  behavior  is  that  as  the  algorithm  proceeds  the  fact  that  the  error  is 
being  measured  in  an  asymmetric  norm  causes  problems.  The  user  of  the  algorithm  should  be 
careful  to  note  that  while  in  principle  all  norms  on  finite  dimensional  spaces  are  equivalent  the  use 
of  a  symmetric  norm  results  in  much  better  observed  behavior  of  this  particular  algorithm. 

To  demonstrate  this  hypothesis  the  problem 


LhU  =  F 


is  transformed  into  the  equivalent  symmetric  problem 

D  1  LhD{irxU)  D  l  F  =--  F. 

The  matrix  D  is  a  positive  diagonal  matrix  whose  entries  are  given  by 


(13) 


1,  d. 


a 


-d, 


lx  I 


Applying  this  transformation  to  /,/,  results  in  l/h  which  is  now  symmetric.  The  matrix  L’h  is 


denoted 


Li 


h  ,  d,  ]/</,  d,  | /(/, 


Denote  by  algorithm  A fC!f  algorithm  M(I  applied  to  problem  (13).  1  he  error  resulting  from 

applying  algorithm  M(!“  is  related  to  the  error  from  applying  algorithm  A f(r  as  follows: 
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Lemma  5.1  If  en  is  the  error  after  the  nth  iteration  of  applying  algorithm  MG  then  the  new  error, 
£n+1,  is  related  to  the  error  in  the  symmetric  problem  by 


(. MGs)(D~le )  =  D~l  MGe  =  D~lcn+X 


Proof:  By  the  definition  of  tn , 


tn  =  U  -  Un 


=  U  -Ur 


MGcn. 


D-yen  =  D-XU  -  D~lUn. 

Applying  each  step  of  algorithm  MG*  to  D~1Un  results  in  D~xUn+l  as  a  straightforward  calcula¬ 
tion  shows.  Hence 

MG*  D~xUn  =  D-'MG  Un. 

Applying  this  transformation  to  the  error  cn  and  computing  the  rates 


D~lRn 


results  in  the  more  usual  behavior  displayed  by  the  dashed  line  in  Figure  1.  Figure  2  displays  the 
location  where  the  maximum  norm  of  the  error  is  taken  on  versus  the  iteration  number.  Notice 
how  as  the  iteration  proceeds  the  location  where  the  maximum  occurs  ‘drifts.’  Applying  D”1  and 
then  measuring  lln  has  the  effect  of ‘fixing’  the  location  where  the  maximum  error  is  taken  on. 
Unfortunately  the  entries  of  the  transformation  matrix  D  satisfy 

Dx—*  oo 

as  e  — *  0  so  this  is  not  a  stable  method  for  solving  this  problem  in  general. 

In  summary,  for  problem  (2)  where  the  coefficient  matrix  was  not  poorly  conditioned  but  was 
non-symmetric  the  only  computational  difficulty  was  in  choosing  the  norm  in  which  to  measure 
the  rate  of  improvement  of  the  algorithm. 

5.3  Turning  point  Problem  -  1 

In  this  section  consider  the  solution  of 

euc"  +  (x  ~)u/  0,  (M) 
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Table  4:  Condition  Number  Problem  TP-1 


with 

«e(0)-l,  «t(l)  3. 

This  problem  has  two  boundary  layers,  at.  x  0  and  at  x  -  1.  The  asymptotic  solution  satisfies 

(see  'KreiT-la  and  DorrlOa'  ) 

lirn  u£(x)  2. 

t  *o 

For  the  discrete  problem  again 

LlU  -  F  (15) 

is  solved  where  l/h  is  given  by  equations  (5)  and  (6).  The  condition  number  of  lTh  tends  to  infinity  as 
e  -  »  0  for  a  fixed  h.  'fable  1  displays  the  L1NPACK  estimate  RCOND  for  lTh ,  h  1/64  and  e  —  1 , 
.01,  .005,  .00,’i,  and  .001.  It.  is  important  to  note  that  although  the  operator  corresponding  to  the 
original  discretization  on  the  h  grid  is  ill  conditioned,  h  3,4,...  are  not  ill  conditioned.  In 
fact  in  the  limit  ing  case  where  t  he  coarse  grid  has  one  point  the  system  of  equat  ions  to  be  solved  is 
a  1  x  1  system  which  has  condition  number  1.  This  feature  of  the  multigrid  algorithm  is  reassuring 
to  the  user  since  it  means  that  the  actual  system  being  solved  directly  in  step  1  of  algorithm  MCI 
is  well  conditioned  and  no  special  measures  need  to  be  taken. 

Although  the  coarse  grids  used  in  the  algorithm  do  not.  resolve  the  boundary  layers  the  algorithm 
still  converges.  This  is  because  the  role  of  the  coarse  grids  is  to  solve  for  the  error  in  the  solution, 
rather  than  to  represent  the  solut  ion  itself.  The  boundary  layers  are  not,  seen  in  t  he  expression  for 
the  error.  The  observed  rate  of  convergence  was  independent  of  the  value  of  e  us«'d 

For  all  tested  values  of  e  the  multigrid  algorithm  converged  to  the  solution  computed  by  Gaus- 
sian  Elimination  to  the  algebraic  problem  However,  the  ill  conditioning  of  the  algebraic  system 
resulting  from  the  discretization  of  the  continuous  problem  opens  up  the  question  of  to  what  so¬ 
lution  does  the  multigrid  algorithm  converge?  Since  the  original  system  is  ill  conditioned  there  is 
a  family  of  solutions  ( 1  for  which  the  residual  is  small  Indeed  for  t  .001  the  condition  number 
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Table  5:  Condition  Number  -  Problem  TP-2 

estimate  is  .71  ■  10  *'  which  is  less  than  the  unit  roundofT  of  the  VAX-11/780.  Thus  there  is  rio 
reason  to  expect,  that  the  solution  calculated  by  solving  equation  (15)  is  an  adequate  approximation 
to  the  solution  of  the  continuous  problem.  In  the  experimental  runs  that  was  indeed  the  case  and 
both  algorithm  M(!  and  (Prussian  (elimination  returned  0  for  the  solution. 

Since  for  this;  problem  the  error  is  not  skewed  as  in  problem  (///.),  it  is  not  surprising  to  see 
that  the  ratio  of  improvements  Pn  are  bounded  above  by  the  predicted  rate  from  Table  1  This 
can  be  seen  in  figure  3,  for  e  .001. 

5.4  Turning  point  Problem  2 

In  this  section  the  solution  to 

suf"  (x  0,  (10) 

is  discussed.  '1  he  asymptotic  solution,  see  Krei74a  is 

1 

uf(j-)  1 ,  0  •  x  •  - 

uf(x)  ‘  •  r-  1 

Although  one  sided  differencing,  in  this  case  /,/ ,  is  used  to  compute  the  solution,  the  system  of 
equations  that  is  solved  is  not  ill  conditioned.  As  m  Section  5.3  the  UNPACK  estimate  RCONT) 
of  the  condition  number  was  computed,  see  Tabic  5 

It,  appears  that  the  ( ondit  mn  number  is  related  to  the  boundary  data  For  problem  (  TP- 1 )  t  he 
discretization  Uses  information  from  the  boundary  layer  to  estimate  ue(x )  at  interior  points.  The 
boundary  data  for  problem  {TP-2)  is  well  represented  (no  boundary  layers)  m  the  interior. 

From  the  standpoint  of  the  multigrid  algorithm  again  as  the  grids  get  coarser  the  condition 
number  improves.  Also,  as  for  problem  (TP-1)  the  rate  of  convergence  of  the  algorithm  was 


rgence 


independent  of  the  value  of  e,  and  was  identical  to  the  rate  predicted  in  Section  3.1.  In  addition 
Rn  was  bounded  above  by  the  predicted  rate. 


6  Comparison  of  ljLh  L  p  I ^  and  L-ik 

From  a  computational  point  of  view  one  would  like  to  choose  for  the  coarse  grid  operator  L^h  the 
tridiagonal  operator  analogous  to  Lh  obtained  by  finite  differences.  This  choice  for  L2/1  will  be 
denoted  lI^.  For  now  only  the  symmetric  problem 

-  (pu')'  =  /,  p(z)  >  Po  >  0  (17) 


with  boundary  conditions 

u(0)  =  u(l)  =  0 

is  considered.  Here  p(x)  is  a  smooth  function. 

Lemma  6.1  The  operator  [J,^  is  related,  to  Lih  by 

Uh=  L{dh  +  h2\  -A,  Bj  -C}  } 
where  A]r  C}  depend  on  p(x),  p'(x)  and  p"[x),  and  B}  =  A}  +  Cj. 
iof:  Recall  that  the  coefficients  of  L^h  are  given  in  equations  (10-12).  In  particular 


a  k  -  ~ 
*  2 


Q2/ca2fc-l 

!hk- 1 


Expanding 


<»2  *  ^  P(~r~h), 

(*-:i t  1  p( — ~-h)  and 

,  ,  dk  1  dk  3  . 

thk  1  /»(— — h)  +  p( — - — «) 


in  Taylor  series  about  the  point  x  -  (dk  -  l)/i  and  collecting  terms  yields 
1  <*2k<*-:k  1 


2  A 


ink  - 1 


1  [p(x)  +  |p'(x)  +  jp"(i)  +  Ci/i3][p(x)  -  |p'(x)  +  jp"(i)  -  c2h3} 


2  (p(^)  +  | p'(*)  +  t-p"^)  +  +  (p(i)  ~  £ p'M  +  fp"(x)  “ 

Ip^  +  ^jl^^p^-lCp^l  +  DA3 

2  2p(x)  +  £p"(x)  +  £h3 


=  ~  P{x)  +  h2Ak 

4 


—  cx{  +  h2  Ak 


where  Ak  depends  on  p(x),  p'(x)  and  p"[x). 

A  similar  calculation  yields 

7*  =  7  ld  +  h2Ck. 


In  addition 


otk  +  7  k 


—  c*k  +  h2Ak  +  7/t  +  h  Ck 


=  P!kd  +  h2Bk 


The  error  analysis  in  Section  3.1  shows  that  if  L2k  lS  not  used  then  the  coarse  grid  correction 
will  not  equal  w2k.  Thus  a  bound  on  the  eigenvalues  of  the  eigenvalue  problem 


A  O  =  l2hi’ 


is  needed. 


Lemma  6.2  The  eigenvalues  X  of 


X  L[dhxf=U^ 


satisfy 


. ,  ,  max  1C,] ,  2 

A  -  1  <  —  -fh2  <  oo. 

rrun  |'y,| 


Proof:  By  Lemma  6.1 


Lih  ~  f'2 h  +  ^2|  ^ Cj  l> 


(A  1  )(Lf2dhtP,rP)  =  h2([  -A,  B,  -Ci  ]0,0). 


AY7-. 

>  V 

tV. 


.  .V 

v' 


'Von 

/  VV  i 

.  S  V  N 

>  v 


•.  V  ^  X.  . 


Summation  by  parts  yields 


\  i  .id  ~A'  B *  ~C<  ]V>,  V») 

A  -  1  =  n  - r - 

h2  Ec,(^  -  t/’.+i)2 

-  v,.+i)2 

<  I 

min| 7/  | 

=  tf/12 

<  00 

because  of  the  smoothness  of  p  and  since  \^{d\  >  po  >  0.  In  addition,  since  is  positive  definite 

(A  -  1  ){Lf2fy,rP)  >  0 


which  implies 


Therefore 


A  -  I  >  -  Kh2. 
| A  —  1 1  <  Kh2. 


Recall  from  equation  (9) 

and  if  L^h  is  used  then 


In  the  case  where  L ^  is  used  to  solve 


c°  -  r?  +  l£hU>2h 


(  -  T). 


h[d^>2h  -  L2hU>2h 


the  error  is  denoted  t  .  Now 

il  =  p  +  /2\(tu2h  -  ip2h)- 

Since  Range  l£h  and  Nullspace  ij^Lh  are  Lh — orthogonal  (see  [Kamo85a]  ), 


Applying  Lemmas  6.1  and  6.2 


\\I2h{w2h  ~  02/i)|iih  =  (Lhl2h(w2h  ~  4>2h),  l2hiw2h  ~  'l>2h)) 

-  \(L2h^2h  ~  ^2h),(W2h  -  ^2h)} 

-  ^(L2h[u>2h  ~  {L^i)  1  L2KW2h),{w2h  -  {L[i)  L2hW2h)) 
=  \{FLl£w2h,FLX£w2h) 


where 


The  eigenvalues  of  L2^( L 2^)  are  the  same  as  the  eigenvalues  of  ( L 2^)  L2h ,  so  by  Lemma  6.2 


\  (FLl2'h2w2h,FL\/h2w2h)  <  -(Kh^CLtfw^LtfrviH) 


=  )  ( L2hU>2h,  tn;/)) 


=  {Kh2)2(LhI^w2h,[^w2h). 


Therefore 


-  V'2a)|ltfc  -  (^^2)  II ^2/1^221  II 

l|i‘llL  <  IMlL  ^  (^^)2l|/2W!lL- 


and 


since  ll^i!^  <  1 


Finally,  the  rate  using  Z,^  is 


l-u 

!■>. 


i  Kfr. 


Recall  the  rate  using  L2h  is 


Remark:  Although  the  discussion  in  the  previous  section  is  restricted  to  symmetric  problems,  the 
same  result  extends  to  the  general  two  point  boundary  value  problem 

-  ( pu')'  t-  bu'  +  cu  =?  /,  p(x)  >  po  >  0,c(x)  >0  (18) 


since  problem  (18)  can  be  rewritten 


-  (pu1)'  t  cu  =  / 


where 


The  function  g(x)  is 


P 

c 


f 


9(z) 


c[x)_ 

p(x) 


!?(x) 


M 

p{x) 


9[x) 


b{x)  -  p'(x) 

p(x) 


dx). 


Note  that  p,  r  and  /  are  defined  since 


(19) 


(20) 

(21) 

(22) 


p(x)  >  p(,  >  0. 


G.l  An  Illustration 

As  an  illustration  of  the  effect  of  replacing  L^h  with  consider  using  for  the  smoother  C  in 
algorithm  Md  one  sweep  of  odd  (lauss-Seidel  smoothing.  In  other  words 

(d  r)t  l\  for  i  0  mod  2 

and  for  i  1  mod  2 

(<•  f),  i  *  iJ'n  l  *  /• . 


This  guar  ant  (‘os  that  the  error  (  lies  completely  in  Range  I ^  In  other  words,  the  r?  term  is  0  since 


.  |Ui»-  «».  !'*  !*• 


JULlVj 


Alt  *»-J  *-M  4.1 


0  —  {Lhf),  l2hw2h) 

-  {V,  ['hl2hw2k)- 


In  addition  for  points  xt ,  with  i  =  1  mod  2, 


0-  r,  -  Lht  --  ( Lhn)i  (Lh^hV’lh), 


which  implies  r/,  -  0  since  r/t  j  t  j  0. 

Thus 

f  —  l2hw2h 

for  some  If  Z-2h  <s  used  then  one  iteration  of  the  algorithm  results  in  completely  solving  the 

problem. 

For  a  numerical  example  problem  (17)  is  solved  with 

p(z)  eco,(™'\ 

The  right,  hand  side  /  is  chosen  so  that  the  the  true  solution  is 

Utrue(x)  ---■  sill(trx). 

As  expected  when  L->h  >s  used  the  algorithm  converges  in  one  step.  When  is  used  then 
the  observed  rate  of  convergence  corresponds  to  the  Kh 2  term  of  Section  6.  Table  6  displays  the 
observed  rate  of  convergence  versus  h  for  four  different  choices  for  the  number  of  grids.  The  column 
p(h)  corresponds  to  the  observed  rate  of  convergence  for  each  value  of  h.  The  value  of  fi(h)  is 


Since  the  error  in  the  two  grid  case  is  0(h~)  one  expects  that  f'(h)  >  1  as  h  *  0  for  two  grids. 

Indeed  this  is  the  case.  The  reason  why  p(/i)  varies  with  the  number  of  grids  used  is  that  there  are 
‘pollution’  effects  from  not  solving  each  coarse  grid  equation  exactly  In  other  words 

>'  '//.  *  I-ihV2h  t  ...  t  o(tr), 

where,  the  O(h')  term  corresponds  to  the  error  made  by  not  using  L^h- 


'A'-’.Vj'Sv.x  ‘v'v  . 


•   -  -  •  i  % 


/  / 
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2  grids 

3  grids 

4  grids 

r  grids 

h 

p(h) 

6(h) 

P[h) 

6(h) 

pW 

6(h) 

P(h) 

6(h) 

1/16 

.50  x  10  2 

.27  x  10”1 

.85  x  10"1 

1/32 

.12  x  10  2 

4.17 

“62  x  10“2 

4.35 

.28  x  lO”1 

3.04 

.86  x  10-1 

1/64 

.30  x  10  3 

4.00 

.15  x  10”2 

4.13 

.65  x  10-* 

4.31 

.29  x  10" 1 

3.01 

1/127 

.75  x  10  4 

4.00 

.38  x  10  ' 3 

3.95 

.16  x  IQ"2 

4.06 

.66  x  10"2 

4.33 

Table  6:  Rate  using  Ljf 
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by  a  multigrid  algorithm  is  considered.  Theoretical  and  experimental  results  for  a  number  of 
different  discretizations  are  pieseuted.  The  theoretical  and  observed  rates  agree  with  the  results 
devel  oped  iii  an  earner  work  of  Kamowitz  and  Harter. 

In  addition,  the  rate  of  convergence  of  the  algorithm  when  the  coarse  grid  operator  is  the 
natural  finite  difference  analogue  of  the  fine  grid  operator  is  presented.  This  is  in  contrast  to 
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